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ABSTRACT 

An algorithm for creating synthetic telescope images of Smoothed Particle Hydrodynamics 
(SPH) density fields is presented, which utilises the adaptive nature of the SPH formalism 
in full. The imaging process uses Monte Carlo Radiative Transfer (MCRT) methods to model 
the scattering and absorption of photon packets in the density field, which then exit the system 
and are captured on a pixelated image plane, creating a 2D image (or a 3D datacube, if the 
photons are also binned by their wavelength). The algorithm is implemented on the density 
field directly: no gridding of the field is required, allowing the density field to be described to 
an identical level of accuracy as the simulations that generated it. 

Some applications of the method to star and planet formation simulations are presented 
to illustrate the advantages of this new technique, and suggestions as to how this framework 
could support a Radiative Equilibrium algorithm are also given as an indication for future 
development. 
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1 INTRODUCTION 

The physics of dust is of paramount importance when modelling 
star and planet formation: dust signatures are seen in both molec- 
ular clouds and circumstellar discs, and are a crucial component 
in their observed properties over a wide range of wavelengths. In- 
deed, the existence of dust is essential if planets are to be formed 
inside these discs. The presence of dust has important effects on the 
radiative signatures that can be detected by astronomers: dust will 
scatter and polarise at short wavelengths, as well as reprocessing 
this radiation to longer wavelengths. The nature of the scattering 
and polarisation will depend on the geometry of both the system 
and the dust in the system, as well as the physical properties of the 
dust itself (composition, morphology, grain size, etc). 

In recent years, telescopes/networks such as HST, SCUBA, 
MERLIN and Spitzer have efficiently probed stellar systems from 
the IR to the radio, with future ground and space-based mis- 
sions such as Herschel, ALMA, SCUBA II, JWST, SPICA and e- 
MERLIN improving the quality of this data. Astronomers will soon 
be faced with a wealth of new, high-fidelity astronomical data on 
stellar objects across a wide spectral range, allowing detailed stud- 
ies of the evolution of dusty disc systems, in particular the interplay 
between disc dust and disc gas. For numerical simulations to inform 
these observations, the simulation of radiative transfer (RT) in these 
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systems must be pursued, so that imaging of numerical simulations 
can provide theoretical insights to observations. 

Monte Carlo Radiative Transfer (MCRT) has become a 
popular tool in simulating and imaging astrophysical system s 
dWhitnev & Hartmand[l992l ; IWood et alj|l996l : iPinte et alj|200r3) . 
The technique involves the tracking of photon packets in the 
medium, allowing them to scatter and absorb in the medium (as 
well as attaining a non-zero polarisation). The photons are tracked 
until they escape the medium, and can then be captured on an im- 
age plane. Traditionally, the method requires the density field to be 
defined to as high an accuracy as possible - this is usually achieved 
by gridding the density field in 3 dimensions. 

Smoothed Particle Hydrodynamics (SPH) is a Lagra ngian 
method which represents a fluid b y a particle distribution JLucvl 
Il977l : iGingold & Monaghan] 1 19771) . Each particle is assigned a 
mass, position, internal energy and velocity. From this, state vari- 
ables such as density can then be calculated at any posi t ion in 
the s ystem by interpolation - see reviews by iMonaehanl {HH 
2005). It has been successful in modelling astrophysical sys- 
tems of various scales and geometries, from protostellar sys- 
tems to galactic systems and beyond. Its key advantage is the 
ability to follow the change of density through many orders of 
magnitude adaptively. Gravity can be included in SPH calcula- 
tions, and optimised us ing traditional hierarchical tree methods 
dHernquist & Katzlll989l) . The formalism can also be modified to 



2 Duncan Forgan and Ken Rice 



simulate magneto-hydrodynamics (MHD) l Hosking & Whitworthl 
l2004|Price & Monaghanlf2004l2005l ; |Price & Batel2007l) . 

Radiative physics is crucial to the simulation of any astro- 
physical fluid. Radiative transfer in SPH has a long history, with 
efforts rangin g from simple par ametrisations for the radiative cool- 
ing time (e.g . iRic e et al. 20031 ), opti cal-depth dependent radiative 
cooling (e.g. Stamat ellos et al.ll2007l). through to flux-limited dif- 
fusion models (e.g. Whitehouse & Bate 2004; Bastien et al. 2004; 
Viau et all 120061; IWhitehouse & Bate! I2OO6I; Waver et all boOl ; 



Forga n et al. l2009h . While these algorithms are extremely well- 



suited to calculating gas temperatures during runtime, they lack the 
precision and insight offered by radiative transfer techniques which 
do not average over frequency, such as MCRT. 

MCRT techniques can be applied to the SPH density 
field, and have been on several occasions. Early attempts be- 
gan by binning the particle distribution onto a grid - how- 
ever, the choice of geometry strongly influences the final 
gridded field, and often adaptive meshes are required to 
corre ct ly represent the matte r distribution dOxlev & Woolfsonl 
20031; iKurosawa etalJ |2004| ; IStamatellos & WhitwortrJ 120051 ; 



Bisbas et alj|2009l). Later efforts have utilised ray-tracing directly 



in SPH field s 1 Kesse l-Devnet & Burkertl 120001; lAltav et all 120081 ; 
IPawlik & Schavell2008l) . allowing the full power of the SPH for- 
malism to be applied to the calculation of optical depths. However, 
these methods currently only calculate optical depths along the ray 
for the purpose of photoionisation, etc, and do not account for de- 
tailed scattering and polarisation, which are important in imaging 
small-scale systems such as protoplanetary discs. 

This paper introduces an algorithm for imaging SPH simu- 
lations of star and planet formation directly using MCRT, with- 
out requiring gridding, and including scattering and polarisation. It 
utilises the same techniques that a traditional gridded MCRT code 
uses, with the advantage that it can trace rays in the density field 
with the same adaptive capability as SPH, such that it can model 
radiative effects with at least the same resolving power. The paper 
is organised as follows: section [2] will outline the algorithm, sec- 
tions [3] and 4 will describe some applications of the technique to 
imaging numerical simulations, and section [5] will summarise the 
results of this work. 



2 METHOD 

2.1 Monte Carlo Radiative Transfer 

For completeness, a summary of MCRT is given here. Photon pack- 
ets are emitted from sources (these can be point sources or diffuse 
emission from the density field itself). These packets (hereafter re- 
ferred to simply as photons) then travel through the medium, and 
interact with the medium stochastically, reproducing the statistical 
scattering properties of the medium. The exact formalism is sum- 
marised below. 



ulations used in this work produce temperatures for each SPH par- 
ticle self-consistently, the system is assumed to already be in tem- 
perature equilibrium. It is also assumed that the dust is thermally 
coupled to the gas, and that T dus t = T rad i a uon- This is suitable 
for most purposes, except where significant stellar irradiation dom- 
inates the radiation field (not modelled in the simulations described 
in this paper). However, as has been already said, radiative equilib- 
rium techniques can be used to sidestep this issue: an implemen- 
tation of radiative equilibrium in SPH fields is discussed in a later 
section. 

It is reasonable to (initially) assume all objects emit according 
to a blackbody spectrum, which then gives: 



Lstar = 47T 2 i?j 



Lgas — ^Igas 



B„(T s )dv 



e„(r)du 



dr, 



(1) 



(2) 



for source and diffuse emission respectively (given the frequency 
range [v m i n ,v max \ of interest), where e„(r) is the emissivity in- 
terpolated from the SPH particle field, i.e. 



e„(r) = 4tt\ W(r - r h h), 

i Pi 



(3) 



where the sum over j indicates a sum over nearest neighbours, k v is 
the absorptive opacity, Tj is the temperature of each SPH particle, 
rrij and pj are the masses and densities respectively, and W is 
the smoothing kernel (see section l2".2.1l for more detail on the SPH 
formalism). 

There is an extra subtlety regarding the frequency distribution 
of photons emitted from the gas. The frequency distribution will 
depend on the local emissivity, which as shown above is an inter- 
polative sum. While in theory the full form of equation (O should 
be used to calculate gas luminosity (and the frequency distribution 
of the photons it emits), this paper approximates the sum using the 
contribution from one particle only, i.e. 



B v {T z )dv. 



(4) 



This saves computational expense, and is sufficiently accurate for 
the examples in this paper, bearing in mind a) that the approxima- 
tion breaks down only in the inner regions of the disc, which are 
already under-resolved for other reasons (see section |4~2t and b) 
the effective resolution of the images is too low for this effect to be 
significant. 

The luminosity of each object (whether pointmass or SPH par- 
ticle) defines how many photon packets are emitted from that object 
using 



N. 



*f ,objec 



-'object 
Ltot 



(5) 



2.1.1 Emission of Photons 

Photons can be emitted either from a central source (e.g. a pro- 
tostar), or as diffuse emission from the density field itself, which 
is only done if the temperature structure of the diffuse component 
is known (e.g. as outpu t from a simulation). Rad iative equilibrium 
methods in MCRT (e.g. Bj orkm an & W ootfeOQll) can calculate the 
temperature structure directly from any luminosity source (e.g. stel- 
lar, accretion, external radiation fields). However, as the SPH sim- 



2.1.2 Absorption/Scattering of Photons 

In order to correctly reproduce the scattering of photons in the 
medium, the cumulative distribution function (CDF) of interaction 

F(t) = 1 - e~ T (6) 

must be reproduced. In essence, this implies sampling an optical 
depth using 

Tscatter = -1 - log(l - C) ( 7 ) 
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where £ is a random number between and 1. The photon will 
travel a distance L along a constant vector path (hereafter ray), such 
that 



Jo 



(8) 



(where \ v is the total opacity), and then be either absorbed or scat- 
tered. The majority of computation required in any MCRT code is 
the solution of equation (8). As the density field of an astrophysical 
system is in general non-trivial, analytic solutions cannot be ap- 
plied, and numerical procedures must be used. The algorithm must 
be able to calculate the optical depth along the ray for arbitrary 
distances - this requirement is what demands the most CPU time. 
To simplify matters, most MCRT codes use gridding as a means of 
defining the density field simply in space. 

Once the photon has reached the scattering location, its direc- 
tion after scattering must be calculated. The angle of scattering is 
defined by the scattering matrix M which acts on the Stokes vec- 
tor S — (I, Q, U, V), where / is the intensity, (Q, U) are the linear 
polarisations at 45° to each other, and V is the circular polarisation: 



S' = R(tt - i 2 ) M R(-h) S. 



(9) 



The R matrices are Mueller matrices, which describe rotations to 
and from the observer's frame. They are defined as: 



R(4>) 



10 

cos 2ip sin 2tjj 

-sin2V> cos 2ip 

1 



(10) 



The scattering matrix M is dependent on the dominant source of 
scattering in the medium. It can be expressed as a function of sev- 
eral scattering parameters, Mc 



M(6) = a 



Mi 


M> 








M 2 


Mi 














M 3 


-Mi 








Mi 


M 3 



(11) 



The scattering angle O and azimuthal angle <f) are sampled ran- 
domly from the scattering matrix, using the cumulative distribution 
functions: 



F(Q) = 



Where 



J o e Misine'de' 
j; Mi sin Q'dQ' 

Mi - M 2 



Mi + M 2 



P 

— sin 2<f) 



y/Q 2 + U 2 



(12) 



(13) 



(14) 



In general, whether the photon is absorbed or scattered, the total 
number of photons is conserved by forcing absorbed photons to 
be immediately re-emitted. As MCRT deals with photon packets, 
energy can be conserved by reducing the energy of each photon 
packet after a scattering/absorption event by a factor equal to the 
local albedo. This equates with the concept of a fraction of the pho- 
tons in the packet being absorbed by the medium, and the comple- 
mentary fraction being scattered. Other methods can also be used, 
e.g. "killing" a photon if the local albedo is less than a randomly 
sampled value, which saves computing emission from low-intensity 
packets. However, this work uses the former method rather than the 
latter. In the following work, the absorbed photons do not affect the 




Image Plane 



Figure 1. Defining an image plane . 



local temperature structure - it is assumed that the equilibrium tem- 
perature is known (having been produced by the SPH simulation). 



2.1.3 Imaging 

When photons exit the medium, they are captured on an image 
plane, oriented at user-specified angles 8 V , (j> v to the system, at a 
fixed distance d (see FigureQJ. They are then binned by their (x, y) 
coordinates on this plane to provide a pixelated image, averaged 
over solid angle (analogous to imaging in a CCD). If spectra are 
of interest to the user, then these can also be obtained by binning 
in A (or indeed, an entire datacube can be obtained by binning in 
all three). "Classic" MCRT methods do not specify a single viewing 
angle, and instead bin the photons over several lines of sight, which 
allow the construction of a series of image planes from one simu- 
lation. However, such "multi-plane" simulations will have to run 
many more photons than a "one-plane" simulation, to maintain a 
comparably low level of random error associated with each pixeQ 
For this reason the code used in this work adopts the "one-plane" 
imaging scheme. 

This numerical apparatus requires the code to be run once for 
every desired orientation of the image plane. It is possible to avoid 
this issue by defining a series of image planes around the system, 
but the numerical resolution of the Monte Carlo method is reduced 
as the pixel to photon ratio is increased, so this is in general avoided 
to prevent significant computational expense. 

The image plane size is defined in physical units, e.g. an im- 
age that captures features a distance of r au from the centre of 
the system. This prevents photons from features outside r (or pho- 
tons which escape with a vector which does not intersect the plane) 
from being recorded in the final image. Having specified r and the 
distance d, this sets the angular size of the image. Coupled with 
information regarding the angular resolution of the telescope be- 
ing simulated, this defines how well resolved the image is, and the 
calculated flux. In lieu of more sophisticated telescope simulation, 
the image is simply smoothed pixel by pixel with a Gaussian with 
width equal to the angular resolution. While obviously far short of 
a true synthetic image, it remains a useful first-order approximation 
for this work. 



1 This is also a factor when comparing one-plane simulations with different 
pixel resolutions. 
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2.2 Implementing Gridless MCRT 

2.2.1 Optical Depths in an SPH density field 

As has been said previously, the key component of any MCRT code 
(and the source of the greatest computational burden) is the calcu- 
lation of optical depths. This calculation requires a specification of 
the density field at all locations, which can be given by the SPH for- 
malism. The definition of a scalar field A(r) using an SPH particle 
distribution is 

A{T) = ^AMyf( T _ T h) (15) 

i Pj 

where Aj is the value of A at the location of particle j, rrij is the 
mass of particle j, pj is the density of particle j, and W is the 
smoothing kernel (or interpolating kernel). The summation is typi- 
cally over N neig h nearest neighbour particles to the location r. The 
smoothing length of particle j, hj, is defined such that a sphere of 
radius 2hj will contain the N ne i g h nearest neighbour particles to 
j. For example, to calculate density, substitute A — p: 

(0 (r)=J]m i W(r-r fj / l ). (16) 

j 

The sphere that contains the N ne i g h nearest neighbours (i.e. a 
sphere of radius 2hj) is referred to as the smoothing volume. There 
is a subtlety to equation l |16t that should be noted, relating to which 
value of h to use. There are now two means by which to esti- 
mate density: the first is the so-called "gather" method, where the 
smoothing length h — hi is defined for the location r;, and 

p(r i )=Y,m j yf(r i -r j ,h i ) (17) 

3 

Where the index j refers to all particles which are contained within 
a radius 2hi of the location r». 




Figure 2. Illustrating the "scatter" method. Particles only contribute to the 
density along the ray if their smoothing volume intersects it. 



S 




The second meth od (which is used in this work and in SPHRAY 
lAltav et ai]|2008l) is the "scatter" method. The smoothing length 
h — hj is used - in this formalism, the density at any one loca- 
tion is calculated by adding the contributions from particles whose 
smoothing volume intersects the location: 



(18) 



In the context of ray tracing, the density along the ray is affected 
only by particles with smoothing volumes that intersect it (see 
Figure[2}. By determining which particles intersect the ray, the rest 
of the particle distribution can be ignored for the purposes of cal- 
culating optical depth, reducing computational expense (whereas 
with the gather method, the ensemble of particles contributing to 
the calculation changes significantly with position, and requires 
the inclusion of a larger subset of SPH particles to perform the 
calculation). 

Using the scatter method, the column density E along the ray is 



Jo Jo 3 . =1 

Which can be rearranged to give 



3 = 1 



mjW(ri — Tj, hj)di 



(19) 



(20) 



Figure 3. Optical Depths through a single smoothing volume . 

The integral is now decomposed into TV integrals, where TV is the 
number of particles intersected by the ray. Each integral is defined 
by the impact parameter b (see Figure The calculation itself 
can be performed for a smoothing volume of h = 1, and scaled 
upwards (this is due to the construction of the smoothing kernel). 
The entire optical depth calculation has been decomposed into the 
repetition of a single algorithm for calculating the optical depth 
through a single smoothing volume. This calculation will now be 
expounded. 



2.2.2 The Optical Depth Calculation for A Single Particle 

Consider Figure [3] The ray (with direction vector n) intersects the 
sphere with impact parameter b. If the ray penetrates a distance x 
into the sphere (out of a total possible distance s), then the inte- 
gral can be defined analytically, given the functional form of W. 
defining 



f — r/h 
b = b/h 
x = x/h 



+ b 2 



(21) 
(22) 
(23) 
(24) 
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Figure 4. Schematic of an Axis Aligned Bounding Box (AABB). For a 
given set of occupants of a tree cell (solid line), the AABB (dashed line) 
is set to the minimum dimensions required to completely contain their 
smoothing volumes. 

(and ff = rf/h), the integral in equation d20t becomes: 

f f, 2 W(f')df' x < s/2 

I — { r i . (25) 

\ g W(f')df' + f b f W{r')dr' x > s/2 

The column density through a smoothing volume (for any impact 
parameter and any distance into the sphere) can now be calculated; 
multiplying by an opacity then provide s an optical dept h. Typically, 
W is constructed using cubic splines (M onas hlmll 19921) . The kernel 
used in this work is 

f 1 " ¥ 2 + I? 3 f < 1 
W(f)=i \(2~r) 3 Kf<2 (26) 

{ f > 2 

This provides compact support (i.e. it reduces to zero outside the 
smoothing volume), and is simple to integrate. 

2.2.3 Ray Tracing in an SPH density field 

With a prescription for calculating optical depth for a single sphere 
in place, a scheme for calculating ray/sphere intersections must be 
constructed. To this end, the code creates a data object called a 
raylist, which stores (in order of intersection) all particles that the 
ray (given its origin and direction vector) will intersect. Once the 
list is created, the optical depth can be calculated quickly using 
equation dl8t . 

The construction of the raylist must be computationally effi- 
cient for the cod e to be effective. The procedure is similar to that 
implemented bv lAltav et al. ( 2008) in SPHRAY; the code constructs 
an octree to spatially index the particles efficiently (as there may 
be density changes over several orders of magnitude). The cells 
either contain child cells, or particles (the leaf cells). The tree is 
constrained to have a maximum number of particles in each leaf. 
All cells have an associated Axis Aligned Bounding Box (AABB), 
which is the minimum box size, aligned to the three cartesian axes, 
to contain all the smoothing volumes of the particles in the cell (see 
Figure[4]l- These AABBs are necessary as tree nodes may contain a 
particle, but not its entire smoothing volume 

This allows the determination of intersections between the ray 
and the cells (or more correctly, their AABBs). Starting with the 
root cell, each child cell is tested for intersection, constituting a 
walk through the tree. If a leaf cell is intersected by the ray, then 



the particles in the leaf are tested for intersection (by calculating 
their impact parameters). This ensures that only a minimum frac- 
tion of the particles in the system need testing for intersection. This 
illustrates the necessity of AABBs; tree nodes may contain a par- 
ticle, but not its entire smoothing volume. Thus, calculating inter- 
sections between a ray and tree nodes may miss contributions to 
the density field from smoothing volumes that cross node intersec- 
tions. Tests for intersections betwee n rays and AABBs are carried 
out using the ray slopes algorithm dEisemann et alj|2007l) . which 
has been shown to be faster than other commonly used me thods, 
such as using Plticker coordinates dMahovskv & Wvvillll2004l) . 



2.2.4 Determining the Scattering Location 

An important facet of an MCRT code is the determination of the 
scattering location of the photon. In general, the scattering location 
will occur inside a smoothing volume, and possibly at a location 
where the density depends on the contributions from several parti- 
cles. Therefore, when attempting to determine the scattering loca- 
tion, it is important to define four classes of particle: 

(i) Particles that do not intersect the ray {unlisted) 

(ii) Particles that intersect, but do not contain the location of 
emission (distant-listed) 

(iii) Particles that intersect, and contain the location of emission 
in front of them (front-listed) 

(iv) Particles that intersect, and contain the location of emssion 
behind them (back-listed) 

The classes are illustrated in Figure [5] Particles of class (i) obvi- 
ously do not affect the calculation - particles of class (ii) are ac- 
counted for simply. Particles of classes (iii) and (iv) will have dif- 
fering effects on the optical depth calculation, and will require sep- 
arate treatments. 

The scattering location is determined by iteration: firstly, the 
optical depth is calculated particle by particle using the raylist un- 
til the optical depth exceeds the randomly selected optical depth 
Tscatter at particle k. Then, the optical depth is calculated from the 
beginning of the sphere for particle (k — 1) (ensuring that all poten- 
tial contributors before and after this location are accounted for), it- 
erating over distance until the answer converges on r 3catter . As the 
optical depth always increases with distance, convergence can be 
achieved with simple algorithms and relatively little computation. 
This code uses a recursive bisector algorithm to perform the itera- 
tion. Starting from the path length between the beginning of sphere 
(k — 1) to the end of sphere k, this value is halved recursively until 
the correct optical depth is obtained (to within some tolerance) or 
until the path length reaches a minimum value (defined as a fraction 
of the smallest smoothing length in the simulation). 



3 TESTS AND APPLICATIONS 
3.1 Comparison with Analytic Results 

To confirm the raytracing component of the code was working cor- 
rectly, a simple test case was devised (Figure [6}. Consider a uni- 
form density sphere, with radius R, density po and total opacity 
X- A point source is located a distance D from the centre of the 
sphere. It emits rays at an angle 9 from the vector connecting the 
source and the sphere's centre. The optical depth t(6) therefore has 
an analytic solution: 
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X- 






Figure 5. The four classes of particle in the system: unlisted (top left), distant-listed (top right), front-listed (bottom left) and back-listed (bottom right). The 
"x" denotes the emission location of the photon. 




Figure 6. Schematic of the raytracing experiment. The optical depth of the 
ray can be calculated analytically to compare with the code's output. 



t{6) = 2 PoX ^R 2 -D2sm 2 (8) (27) 

Three SPH snapshots were generated, each containing a uniform 
density sphere (mass 1 Mq, R = 2133 au). Three were generated 
to check convergence: snapshot 1 used 10 5 SPH particles; snap- 
shot 2 used 5 x 10 5 ; snapshot 3 used 10 6 . A point source was 
placed at a distance D = 4000 au, and the optical depth along 
the ray (assuming an opacity of unity) was calculated as a func- 
tion of 9. The numerical results are compared with the analytical 
result (scaled such that the maximum optical depth is 1) in Figure 
[7] The column density is subject to the underlying random noise 
(at a level of around 5%) associated with generating an SPH snap- 
shot (which has not undergone any settling). However, the results 
vary little with increasing particle number, showing that the column 
density has converged even for the relatively low particle number 
oflO 5 . 
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Figure 7. Results of the raytracing experiment. The black solid line indi- 
cates the analytical result; the coloured dashed lines show simulation results 
for increasing particle number. 



3.2 A Low Mass Companion for HL Tau? 

Grea ves et ail J2008I) imaged HL Tau using the Very Large Array 
(VLA) with a resolution of 0.08" at a wavelength of 1.3cm (corre- 
sponding to a spatial resolution of 10 au at HL Tau's distance of ~ 
140 pc). They discovered excess emission at ~ 65 au, which they 
identified as a candidate protoplanet in the earliest stages of forma- 
tion, a possible exa mple of prot oplanetary disc fragmentation form- 
ing a bound object dBossll997l) . To lend weight to their hypothesis, 
they conducted SPH simulations (containing 250,000 particles) of 
an unstable star-disc system with similar parameters to HL Tau, in 
which a clump forms at ~ 75 au, with a similar mass as that de- 
duced for the candidate (s ee Figure[8l > . The s imulation uses an SPH 
code based on the work of lBate et al.l dl995T) . It employs individual 
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particle timesteps, and models hydrod ynamics and grav ity, with ra- 
diative cooling as prescribed by Stamatellos et al. (2007). What can 
a telescope like ALMA be expected to see in HL Tau? Will the can- 
didate protoplanet be resolvable at millimetre wavelen gths? 

The SPH simulation from iGreaves et alj (2008) was there- 
fo re imaged, using t he same dust scattering parameters as used 
in IWoodetalJj200afi The wavelength and resolution parameters 
were set to be representative of ALMA0 the wavelength range 
is [0.01, 0.1] cm, with an effective resolution of 0.01". The star is 
0.5 Mq, and the disc is 0.3 Mq, with an initial surface density pro- 
file of E oc r _1 A dust to gas ratio of 0.005 (corresponding to 50% 
solar metallicity) was assumed throughout. Dust sublimation is pre- 
scribed by enforcing any SPH particle with temperature greater 
than 1600 K to have zero dust mass. Te mperatures were cal culated 
using the equation of state o utlined in Forgan et al. ( 2009) (more 
inform ation can be found in IStamatellos et alj 20071 : Bolev et all 
l2007l) . The star's emission was modelled by a blackbody source 
with radius 2Rq and temperature 2000K. 

Figure[8]shows that the clump, which reaches temperatures of 
around 1500 K at its centre is detectable by its emission, although 
somewhat fainter than the main star-disc system. In fact, it should 
be expected that the clump is fainter still, if dust sublimation is 
more appropriately modelled as opposed to a one temperature cut- 
off. The m=2 spiral mode attached to the clump is also apparent; 
however with a flux of around 0.001 mjy/beam, it is unlikely that 
ALMA will be able to image it, highlighting the need for resolution 
of order ~ 1 au at 100 pc and sufficient sensitivity in the detection 
of disc spiral arms. The peak emission is around 0.5 Jy, which is 
within a factor of two of the observed SCUBA measurements of HL 
TauQ. The remaining discrepancy is most likely due to uncertainty 
in the selected disc and star parameters (as well as the dust data). 
To summarise, the clump is indeed detectable with ALMA, but a 
telescope with better sensitivity is required to detect the spiral arms 
correctly. 



3.3 The Observational Features of Stellar Encounters 

It is expected that in an average stellar cluster, at least one star-disc 
syste m will undergo a close encounter with another star in its life- 
time dClarke & Prmgldll 99 ll ; Morgan & Ricell2009bh . The effect of 
an encounter perturbs the disc, modifying its surface density pro- 
file and triggering enhanced accretion. The effect of the enhanced 
accretion produces a stell ar outburst. This feature was linked to th e 
FU Orionis phenomenon teonnell & Bastie rll992l : |Pfalzned200"8l) . 
but it has been shown that such encounters are too infrequent to re- 
produce the correct ob servational signatures of FU Ori outbursts 
(IForgan & Ricejl20q9bl). 

iForgan & Ricd fe009ah carried out a series of SPH simula- 
tions to study the effects of a close encounter on a protostellar 
disc's dynamics and structure. As a second application of the imag- 
in g techniques discussed here, the reference simulation described 
in Forgan & Ricej d2009ah was imaged at two points in its evolu- 
tion: the initial pre-encounter phase, w here the disc remain s in a 
marginally stable, self-gravitating state (Lodato & Rice 2004]); and 
during periastron, where the disc is heated strongly by the sec- 
ondary's motion through it. 

The primary is 0.5 Mq (modelled by a blackbody source with 



radius 2Rq and temperature 2000K), and the secondary is 0.1 Mq 
(modelled by a blackbody source with radius O.2i?0 and temper- 
ature 1000K). The disc is 0.1 Mq, with an initial surface density 
profile of £ oc r _1 (and dust to gas ratio of 0.01, i.e metallicity 
equal to solar). The imaging parameters are the same as those of 
the HL Tau example above ([0.01, 0.1] cm, with an effective reso- 
lution of 0.01"). 

Disc asymmetry plays an important role in the imaging of this 
system (Figures[9land|10t. In Figure[9] the disc displays a non-zero 
ellipticity due to its spiral structure, with shadowing indicating the 
increased surface density of the disc corresponding to the spirals 
themselves. Again, however, these shadows are only apparent if the 
sensitivity of the telescope is 1 ^Jy/beam, well beyond the ALMA's 
current continuum sensitivity estimates. This ellipticity is enhanced 
during the encounter (Figure [10}, with the semi major axis of the 
ellipse aligned with the orbital vector of the primary and secondary. 
The erasure of the strong spiral structure during the encounter re- 
sults in a detectable tidal arm with flux around 1 mJy/beam. The 
stellar emission is dwarfed by the disc at these wavelengths. The 
inner regions of the disc are not hot enough for dust sublimation to 
open a resolvable gap (the intrinsic inner gap is very much below 
the resolution limit). Again, the resolution of spiral structure in a 
self-gravitating disc will not be possible with ALMA unless higher 
sensitivities are achieved, but it will be possible to detect enhanced 
emission from a tidal spiral wave generated as the result of a stellar 
encounter. 



4 DISCUSSION 
4.1 Runtime Scaling 

As the SPH systems being imaged by this code will be in gen- 
eral disordered and impossible to render analytically, a true run- 
time scaling is not feasible. However, an example scaling assuming 
simple geometry can be calculated. 
In general, the runtime T goes as 



T ~ N^Nsteps 



2 Note that the examples shown here approximate \v = K f 

3 http://www.eso.org/sci/facilities/alma/ 

4 http://www.jach.hawaii.edu/JCMT/continuum/calibration/sens/calibrators.html 



(28) 



Where iV 7 is the number of photons emitted by the code, and 
N s teps represents the computational expense required to track one 
photon from emission to capture. This can hence be written 



T ~ AL < AT > N scatt 



(29) 



Where < N ray > is the mean number of particles intersected 
by any ray, and N aca tt indicates how many times a photon will 
scatter before it exits. Generally, N sca tt ~< T > 2 , and 



^ N ra y Np Pintersect 



(30) 



Where N p is the number of SPH particles, and Pintersect is 
the probability that any one particle is intersected by the ray. This 
can be estimated for simple geometries: assuming the SPH particle 
distribution is a sphere, then 

Pinter sect ~ 1 — C Rna 3 (31) 

Where R is the sphere's radius, n — N p /V ~ N p /R 3 is the 
number density of the sphere, and a s ~ 4ir < h > is the average 
cross-section of the SPH particles. In a homogeneous sphere, 



<h>r 



1\3 

9 j 



R 



(32) 



Combining these results gives 
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Figure 8. Left: Surface density plot of the HL Tau simulation to be imaged. Right: image of the HL Tau simulation, integrated over the wavelength range 
[0.01, 0.1] cm. The green cross indicates the pixel with maximum intensity. 
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Figure 9. Left: Surface density plot of the disc pre-outburst. The secondary is out of frame. Right: image of the simulation, integrated over the wavelength 
range [0.01, 0.1] cm. The green cross indicates the pixel with maximum intensity 







Figure 10. Left: Surface density plot of the disc at periastron. Right: image of the simulation, integrated over the wavelength range [0.01, 0.1] cm. The green 
cross indicates the pixel with maximum intensity 
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(33) 



This runtime scaling shares common features with grid-based 
MCRT codes, in particular the dependence on optical depth and 
number of particles/cells. For illustration, the simulations used in 
this work typically took approximately 100 CPU hours to create 



an image, with N p = 500, 000, N ga mma 
significant optically thick components. 



4.2 Resolution 



10 in systems with 



SPH simulations are typically subject to several resolution con- 
ditions, princip a l amo ng them the Jeans Resolution criterion of 
iBate & Burkertl dl997l) to correctly resolve fragmentation. Imag- 
ing these simulations is no different - it must be demonstrable that 
the local density field is sufficiently resolved to allow a satisfactory 
calculation of the optical depth (more specifically the optical depth 
to the photosphere, as this is where most of the received emission 
will originate). In terms of the smoothing length, h, and the mean 
free path A m f p , this condition is 



h < Xmfp 



pK 



(34) 



As the value of h is inversely proportional to the local number den- 
sity of SPH particles, the resolution of an MCRT image of an SPH 
density field is closely linked to the total number of particles in 
the simulation, an intuitive result. This condition is satisfied in the 
outer disc easily, as the mean free path is of order the system scale 
height, (which is resolved by several smoothing lengths). However, 
within the inner 5-10 au of the disc, two separate issues arise: 

(i) The smoothing length exceeds the local mean free path. This 
is partially compensated by SPH's smoothing prescription - the 
density field is defined continously across all space, and therefore 
affords a kind of sub-h resolution. However, this smoothing erases 
information about fluctuations in the density field below this length 
scale, and this will lead to an overestimation of the escaping flux. 

(ii) The disc is not well resolved several scale heights above the 
midplane, where the photosphere is located. This under-resolving 
will result in a decrease in scattering events at the photosphere, 
affecting the scattered light component of the flux. 

These two issues show that inside ~ 10 au, the flux es- 
caping from the disc will be subject to resolution-based error. 
Users concerned with the inner regions of these discs could 
adopt a particle splitting pres cription to boost the resolution (e.g. 
iKitsionas & Whitworthll 20021) . However, as this work is concerned 
mainly with millimetre emission from the cooler, well-resolved 
outer disc, this emission will be adequately resolved, and the is- 
sues described above will not play a significant role. 



4.3 Radiative Equilibrium: A Proof of Concept 

The work described here has directly utilised the temperature struc- 
tures produced in the SPH simulations, but such outputs are not a 
requirement for imaging them using MCRT. The temperature for 
each fluid element can be calculated using MCRT with so-ca lled 
radiative equilibrium methods (e.g. lBiorkman"& Wood 200l|). In 
essence, these methods involve absorbed photons incrementing the 
local temperature field at the absorption point T(rj) by an amount 
AT(rj), and then being re-emitted. As many photons are passed 



through the system and are absorbed and re-emitted at various lo- 
cations, the temperature structure will relax to the correct value. In 
grid-based methods, the procedure of increasing the local tempera- 
ture is simple - the increment AT is added to the grid cell the pho- 
ton is absorbed in. In SPH fields, the situation is somewhat more 
complex. A scalar field X at a location rj is defined as: 



Pi 



(35) 



(Where we have used the "scatter" interpretation). We can substi- 
tute for the temperature increase AT: 



AT(r i )=E^W(r i 



. hi 



(36) 



The left hand side of equation J36t can be calculated using the tra- 
ditional radiative equilibrium methods, and is known. The desired 
quantities are the ATij, that is the individual increases in temper- 
ature each SPH particle must be assigned. Equation J36t therefore 
acts as a constraint on the ATy, but is insufficient to close the sys- 
tem. A method of closure is presented here as proof that the radia- 
tive equilibrium framework can indeed be incorporated - its appli- 
cation is left for future work. Closure can be achieved by specifying 
the following ansatz: 



ATi, 



AijATginj) 



(37) 



Where Aij is some real number, nj is the separation of particle i 
from the location j and g{rij) is some function that satisfies: 



3(0) = 1, 

g(rij >2hi) = Q 



(38) 
(39) 



i.e. particles with smoothing volumes that do not intersect the loca- 
tion j do not contribute. This will satisfy equation d 3 €>! > provided 



V — 4yp(ry)W(ry,/n) = 1 

V Pi 



(40) 



This has recast the problem from solving for ATy to the prefac- 
tors Aij and function g(rij). This may appear to be a sideways 
step: however, giving a specific example is illuminating. To obtain 
a suitable function g(rij), consider the decay of flux Fo along a 
path of optical depth r: 

F = F e~ T (41) 
Using the Stefan-Boltzmann law, F — aT 4 yields 
a(ATij) 4 = a(AT) 4 e- T *i (42) 
Inspecting the above equation suggests the following form for g: 



g{rij) 



M /4 



(43) 



Assuming that Aij = A for all the normalisation condition 
becomes 



A = 



E,^e-,-/ 4 W(r y ,h,) 



(44) 



The system has now been closed with an ansatz that is physically 
motivated, and reproduces expected behaviour (particles with a re- 
duced optical depth to the location will receive a greater tempera- 
ture boost than others in more optically thick regions). 

Other radiative equilibrium m ethods can also be recast in SPH 
form, such as that of lLucvl dl999l) . which computes grid cell tem- 
peratures by counting the path lengths of any photons that pass 
through it during one iteration, retrieving the temperature structure 
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rapidly within a few iterations. By replacing grid cells with parti- 
cle smoothing volumes, it can be seen that this method is indeed 
amenable to SPH fields also. 



5 CONCLUSIONS 

This paper has outlined a means for applying Monte Carlo Ra- 
diative Transfer (MCRT) techniques directly to Smoothed Particle 
Hydrodynamics (SPH) density fields of star and planet formation. 
In doing so, it gives a means for synthetic telescope images to be 
made, allowing the flexibility of SPH simulations to be retained in 
calculating optical depths, scattering and polarisation. This work 
has shown two applications of the method to SPH simulations of 
protostellar discs: non-trivial features in the imaging of these envi- 
ronments are well-traced, and the resulting images give new oppor- 
tunities for theoretical input to observed discs. It can be shown what 
resolutions and sensitivities are required to observe these features. 
ALMA appears to be able to resolve and image perturbed structures 
like tidal arms in discs undergoing close encounters. But, while it 
has sufficient resolution, it cannot image unperturbed spiral struc- 
ture (or the shadows it casts) without an improvement in sensitivity 
by at least a factor of 100. Its inherently graphical nature allows 
the code to be greatly optimised, whether by standard parallelisa- 
tion methods or by the use of Graphical Processing Units (GPUs) 
to handle the ray intersections. 

While used for imaging in this work, the method is not re- 
stricted to this alone. As the algorithm only modifies how the opti- 
cal depth is calculated (in order to work in SPH fields), it can be ap- 
plied to any circumstances traditional gridded MCRT has been ap- 
plied to. This includes radiative equilibrium simulations where the 
temperature structure is calculated from stellar emission, or moving 
beyond continuum emission to perform line radiative transfer cal- 
culations. In effect this places SPH on a par with grids or meshes 
for MCRT techniques. 
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